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When immiscible wetting and non- wetting fluids move in parallel in a porous medium, an insta- 
bility may occur at sufficiently high capillary numbers so that interfaces between the fluids initially 
held in place by the porous medium are mobilized. A boundary zone containing bubbles of both 
fluids evolve which has a well defined thickness. This zone moves at constant average speed towards 
the non- wetting fluid. A diffusive current of bubbles of non- wetting fluid into the wetting fluid is 
set up. 

PACS numbers: 47.20.Ft,47.56.+r,47.54.-r,89.75.Fb 



When a fluid displaces another one in a porous medium 
the interface separating the two fluids may become un- 
stable. In the case of two-phase immiscible displacement, 
local capillary barriers on pore-scale levels affect the be- 
havior on larger scales, and it turns out that there is 
an extraordinary richness to the ways instabilities occur 
and how the separating interface subsequently develops. 
Depending on several flow properties like viscosity ratio, 
wetting properties with respect to the porous medium 
and how fast the displacement occurs, a wide range of be- 
haviors are found in both drainage and imbibition rang- 
ing from pure invasion percolation to viscous fingering 
p], 0| . A huge effort has gone into classifying and un- 
derstanding this rich behavior, both from a fundamental 
scientific point of view, but also due to its importance in 
a number of very important fields ranging from oil recov- 
ery, to spreading of pollutants in the ground water, to 
problems related to CO2 sequestering. 

It is then surprising to discover that the related prob- 
lem of immiscible fluids flowing in parallel to the inter- 
face between them rather than normal to it in a porous 
medium, has received very little attention in compari- 
son. Such parallel flow is e.g. seen in connection with 
fully developed viscous fingers J3[ and in connection with 
flow in stratified reservoirs 0, U UWM- When the flow 
rate is low so that capillary forces dominate at the in- 
terface, the parallel interface is stable and each phase 
behaves as in a single-phase flow system. Nevertheless, 
it has been recognized that above a certain threshold in 
the flow-rate, but where capillary forces still dominate, 
imbibition processes become important in the evolution 
of the interface and hence the crossflow of the immiscible 
fluids. However, at larger flow rates, where viscous forces 
dominate, shear-driven Kelvin- Helmholtz type instabili- 
ties are believed to occur 0, [TO, Both theoretical 
and experimental work has been invested in studying the 
Kelvin-Helmholtz instability in vertical Hele-Shaw cells 
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[E EH, [HI, as it provides a model for parallel flow in 
porous media in the viscous regime. 

It is the aim of this Letter to investigate the instabil- 
ity that occurs at the interface between two immiscible 
fluids flowing in parallel in a regime where capillary ef- 
fects cannot be ignored. This regime has remained es- 
sentially untouched in the literature. We find that above 
a threshold flow rate and with a viscosity ratio between 
the two fluids favoring the formation of viscous fingers, 
the interface becomes unstable, and a boundary zone ap- 
pears containing intermixed bubbles of both fluids. This 
boundary zone has a well-defined width and moves at 
constant average speed towards the non-wetting fluid. A 
diffusive current of bubbles of non- wetting fluid into the 
wetting fluid is set up, but the situation of bubbles of 
wetting fluid entering the non- wetting fluid is absent. 

This instability may prove to be an important mecha- 
nism for mixing non- wetting fluid into wetting fluid. A 
practical application may be CO2 sequestering in porous 
rock formations. The less wetting gas is blown into the 
porous medium which is already saturated by a more wet- 
ting fluid. A mixing zone will then form at the boundary 
between the gas and the fluid where gas bubbles will be 
generated. These bubbles are then transported into the 
wetting fluid where they eventually are absorbed. 

We study this instability here using a two-dimensional 
network simulator first developed by Aker et al. [14] with 
later extensions by Knudsen et al. [15| and Ramstad and 
Hansen [16[ . The network forms a square lattice oriented 
at 45° with respect to the overall flow direction. Each link 
forms an hour-glass shaped tube. Disorder is introduced 
in the model by having the average tube radii r be drawn 
from a flat distribution on the interval r G (0.1£, 0.4^), 
where £ is the link length. Capillary pressure in the links 
is caused by the presence of interfaces in them. 

As the tubes are hour glass shaped, the capillary pres- 
sure difference caused by a meniscus at position x, the 
distance from one of the two nodes it is attached to, is 
given by p c oc 1 — cos(27nr/-Q. 

We assume cylindrical tubes so that the flow rate q 
in a tube is given by the Hagen-Poiseulle relation from 
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where Ap is the pressure difference between the nodes 
connected by the tube. The effective viscosity is the 
volume-weighted average of the viscosities of the fluids 
contained in the tube. The sum runs over the number of 
meniscii in the tube. We accept up to 10 meniscii in any 
given tube. If this number is exceeded or the distance 
between two bubbles is too small, we merge the mensicii. 

The flow equations are solved by assuming flux con- 
servation at each node, i.e., invoking the Kirchhoff equa- 
tions. This is done by defining a pressure p at each node. 
We use the conjugate gradient method for this [17[. After 
the node pressures have been determined, the positions 
of the meniscii are integrated forwards by an adaptive 
time-step At so that no single meniscus movement ex- 
ceeds one tenth of a tube-length £. When meniscii reach 
the ends of a tube, they are moved into the other eligible 
tubes connected to that node. For details, see [151 ]. 

The flow of the two fluids in the network is controlled 
by the ratio between capillary and viscous forces at the 
pore level and quantified through the capillary number 
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where \i is the largest viscosity of the two immiscible 
fluids, E is the cross-sectional area of the network and 
Qtot is the total flux through this area. 

Besides the capillary number, the ratio between the 
viscosities of the two fluids forms the second dimension- 
less number to control the flow, 



M 



/iw 



(3) 



We set M = 1 in the following. Hence, there is initially 
no pronounced shear in the flow patterns in the network. 

We implement periodic boundary conditions in the av- 
erage flow direction [TBI, [l6[ • This implies that the flow 
configurations sees no boundaries in the flow direction, 
and the fluid configurations may develop over large times 
and distances. There is no periodicity in the direction 
normal to the average flow direction. The boundaries 
parallel to the average flow direction are in contact with 
a reservoir of either wetting or non-wetting fluid. A con- 
stant pressure drop AP is set up across the network in 
the average flow direction, causing the total flux Qtot- 

The network is prepared with either a band of non- 
wetting fluid parallel to the average flow direction, sur- 
rounded by wetting fluid or vice versa. Hence, the sat- 
uration of non-wetting fluid, 5 nw and wetting fluid, S w 
is non-uniform. We show in Fig. [1] the network initially 
prepared with a band of non- wetting fluid in the middle. 
If the pressure drop AP is too small, the interfaces in 
the tubes forming the boundaries between the two fluid 
types will be stabilized by the capillary pressures and the 





Figure 1: Different stages of the development of an initially 
straight band of non- wetting fluid (black) inside a region filled 
with wetting fluid (white). The flow is from top to bottom 
with periodic boundary conditions this direction. The bound- 
aries are open in the transverse direction. The size of the 
network is L x x L y — 64 x 64. 



boundaries are stable. However, when the pressure differ- 
ence is above a minimal value so that the initial capillary 
number Cai n it > Ca m i n , the boundaries destabilize and 
the system evolves. 

Early in the evolution of the system, fingers of non- 
wetting fluid form when the viscosity ratio between the 
two fluids allow this. This is a signature of unstable non- 
wetting front propagation in the viscous regime. The 
fingers are bent in the direction of the average flow. Due 
to the flow typically being at an angle compared to the 
fingers, they are susceptible to break up. The broken off 
fingers form bubbles that migrate into the wetting fluid, 
and consequently the wetting fluid also migrates into the 
non-wetting fluid. This is due to the appearance of an 
effective pressure gradient APj_ normal to the average 
flow direction across the boundary region between the 
two fluids. This gradient is in turn due to the appearance 
of a gradient in the effective permeability. The effective 
pressure gradient APj_ leads to imbibition of the wetting 
fluid into the non-wetting region. This process creates 
a compact front and a saturation profile moving in the 
direction normal to the average flow direction, resembling 
that of Buckley-Leverett flow [H . 

There is a length scale A associated with the saturation 
profile. We define it through the width of the bell-shaped 
non- wetting saturation gradient as shown in Fig. [2] based 
on an average over five samples. From the motion of 
the two maxima of dS nw /dx along the x-axis, which is 
the direction normal to the average flow direction, we 
determine the mean velocity of the non-wetting satura- 
tion profile. The data collapse shown in Fig. [2k, where 
dS nw /dx is plotted against (x — vt)/L x , shows that the 
mean velocity v of the profile is constant and the shape of 
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Figure 2: (Color online) (a) S n w vs. x for different t, (b) 
dSjxw I dx vs. x for different t, (c) dS^/dx vs. (x — vt)/L x 
for different t for a moving front with starting point xo — 0. 
All figures are for an Lx x Ly — 128 x 32 lattice with open 
boundaries in the direction parallel to the overall flow and 
with Cainit = 0.03. 

the profile is also constant. The length scale A is for this 
system X/L x w 0.04. Hence, unlike the Kelvin-Helmholz 
shear instability, the interfacial instability between two 
different fluids in a porous medium proceeds through the 
creation of a well-defined saturation profile characterized 
by a length scale A and an average speed v which remains 
constant. 

The shape of the non-wetting saturation profile corre- 
sponds to a there being a boundary region where bubbles 
are created. This boundary region moves into the non- 
wetting zone. There are no bubbles migrating into this 
zone ahead of the moving boundary region. On the other 
side, there is a diffusion current of bubbles of non- wetting 
fluid into the wetting zone. The diffusing bubbles stem 
from the non- wetting fingers that break off due to the av- 
erage flow being at an angle with respect to the fingers. 

When the two approaching boundary regions eventu- 
ally meet, the middle non- wetting band is destroyed as 
shown in the last picture of the sequence shown in Fig. 

m 

We now reverse the initial configuration so that a band 



of wetting fluid is surrounded by non-wetting fluid. We 
show in Fig.[3]the evolution of the non- wetting saturation 
profile as a function of time. As before, boundary regions 
where bubbles form are created. However, after some ini- 
tial time, they stabilize and do not move. This is in sharp 
contrast to the previous situation where the boundary re- 
gions move with constant mean velocity. This is caused 
by there being no bubble transport outside the bound- 
ary region and into the non-wetting region. Inside the 
wetting band, there is diffusive bubble transport, but as 
the width of the band is finite and it is surrounded by 
bubble-generating boundary regions on both sides, the 
net diffusive current stabilizes at zero. 

We consider in the following the evolution of the total 
flow rate Q tot as the system evolves for both configura- 
tions we have studied. We consider first the case of a 
non-wetting fluid band surrounded by wetting fluid. As 
the flow is sustained by a constant pressure drop across 
the network in the average flow direction, AP, the total 
flow rate Q tot will at all times be proportional to the per- 
meability of the network. We show in Fig.[4]the evolution 
of the total flow rate as a function of time for networks 
initially prepared with a non- wetting band in the middle 
and with a wetting band in the middle. We analyze first 
the case when the networks starts with a non-wetting 
band in the middle. The evolution of Qtot for two dif- 
ferent pressure drops are shown in Fig. [H We see that 
for both pressure drops, Qtot decreases linearly in time 
after an initial transient. In the case of the larger pres- 
sure drop, the red curve, the flow rate starts increasing 
again reaching essentially the flow rate it had initially. 
This behavior can be understood as follows. After the 
initial transient and before the rapid increase of Qtot in 
the high-AP case, the system consists of three zones: 
(1) a non- wetting zone characterized by an effective local 
permeability fc nw ? (2) a boundary zone characterized by a 
local permeability k\ and (3) a zone where non-wetting 
fluid forms diffusing bubbles in the wetting fluid. The 
local permeability here is fc m i x . If the width of the non- 
wetting zone is ^ nw , of the boundary zone is A and of 
the mixed zone is £ m i x , then the total permeability of the 




)i i i i i i i i i 

0.2 0.4 0.6 0.8 



Figure 3: Non- wetting saturation profiles for the initial config- 
uration of a band of wetting fluid surrounded by non- wetting 
fluid. The profile stabilizes at different levels for different 
AP. For higher AP the saturation 5 nw in the boundary re- 
gion stabilizes at a higher level and the wetting front advances 
further. 
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Figure 4: (Color online) Qtot as a function of time for two 
different constant pressure differences AP when the middle 
band is non-wetting. Inset (a) shows the same for one con- 
stant pressure difference when the middle band is wetting. 
Inset (b) shows the flux profile normal to the flow profile for 
an initial non- wetting band in the middle. 



network is given by 

7 — / ^mix ~i 2 A 4w / ,\ 

^eff = ^mix j r K\y r r>'nw ~J • (4 J 
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As the boundary regions moves with constant average 
speed we have that £ m - lx = ^ m ix,o + vt and since L x = 
2^mix + 2 A + £ mw , it follows that £ nw — L x — 2£ miX)0 - 
2A — 2vt. Inserting these two equations in Eq. (|4|) gives 

2vt 

^eff ^eff,0 j [^nw — ^mix] • (5) 
L*x 

Since the viscositities of the two fluids are equal, Qtot ^ 
fc e ff- As k nw is larger than fc m i x , the total flow rate Q to t 
falls off with time. 



Fig.[4]shows that at a larger pressure difference, the to- 
tal flow rate starts increasing again after the linear regime 
we have just described. This is due to the non- wetting 
band in the middle having been depleted and the non- 
wetting bubbles are diffusing out of the network. Then 
the network is being depleted of non-wetting fluid and 
hence interfaces which lowers the effective permeability. 
The opposite situation, a middle wetting band, is shown 
in inset (a) in Fig. [H We see that the total flow rate sat- 
urates, indicating that the system enters a steady state 
as already discussed in connection with Fig. [3l 

We have in our numerical experiments kept the pres- 
sure difference across the network constant. If we rather 
had kept the total flow rate Q to t constant, the instability 
that sets in when Ca > Ca m i n will be much more vio- 
lent. This is so since the pressure drop AP will increase 
to keep Qtot, leading in turn to an acceleration of the 
boundary region. 

In this Letter we have investigated the stability where 
two different fluid flow parallel to each other. We find 
that under constant pressure conditions, for a sufficiently 
high capillary number a boundary region develops with a 
well-defined width when the viscosity ratio between the 
two fluids favor the formation of viscous fingers. This 
region, which essentially is foam, moves at a constant 
average speed into the non-wetting region. On the wet- 
ting side, a diffusive current of non- wetting bubbles away 
from the boundary zone develops. It would be of great 
interest to see this instability reproduced in the labora- 
tory, e.g. in two-dimensional glass-bead filled Hele-Shaw 
cells. 
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